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ABSTRACT 

We perform three-dimensional smoothed particle hydrodynamics simulations in a realistic cosmo- 
logical setting to investigate the expansion, feedback, and chemical enrichment properties of a 200 Af© 
pair-instability supernova in the high-redshift universe. We find that the SN remnant propagates for 
a Hubble time at z ~ 20 to a final mass- weighted mean shock radius of 2.5 kpc (proper), roughly 
half the size of the H ii region, and in this process sweeps up a total gas mass of 2.5 x 10^ M©. 
The morphology of the shock becomes highly anisotropic once it leaves the host halo and encounters 
filaments and neighboring minihalos, while the bulk of the shock propagates into the voids of the 
intergalactic medium. The SN entirely disrupts the host halo and terminates further star formation 
for at least 200 Myr, while in our specific case it exerts positive mechanical feedback on neighboring 
minihalos by shock-compressing their cores. In contrast, we do not observe secondary star formation 
in the dense shell via gravitational fragmentation, due to the previous photoheating by the progenitor 
star. We find that cooling by metal lines is unimportant for the entire evolution of the SN remnant, 
while the metal-enriched, interior bubble expands adiabatically into the cavities created by the shock, 
and ultimately into the voids with a maximum extent similar to the final mass-weighted mean shock 
radius. Finally, we conclude that dark matter halos of at least Afvir ^10^ must be assembled to 
recollect all components of the swept-up gas. 

Subject headings: cosmology: theory — galaxies: formation — galaxies: high-redshift — H ii regions 
— hydrodynamics — intergalactic medium — supernovae: general 



1. INTRODUCTION 

One of the main goals in modern cosmology is to 
understand the formation of the first stars at the end 
of the cosmic dark ages and how they transformed 
the homogeneous primordial univer se into a state 
of ever increasing complexity fe. g.. iBarkana fc Loebl 
2001; Miralda-Escudc 2003; ' Bro mm fc LarsonI l2004t 
Ciardi fc Ferrara ■2005.: .Glover. ,200511 . These so-called 
Population III (Pop HI) stars, predicted to be very 
massive, with M^: > 100 Mq, and forming in DM 
"minihalos" w ith virial masses Myjr > 10^ M © (e.g., 



iKitavama fc Yoshidal[2005h . 

The first SNe exerted important chemical and 
mecha nical feedback e f fects on the early universe 
(e.g., iCiardi fc Ferraral I20050 . First, the charac- 
ter of star formation changed from an early, high- 
mass-dominated (Pop HI) mode to a more nor- 
mal, lower mass (Pop II) mode, once a criti- 
cal level of enrichment had been reached, the so- 
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Bromm et al.1 11 999. 2002; Nakamura fc Umemura 2001 
Abel et all 12002: Yoshida et al.. ,20061: iGao et al. 2007 
O'Shea fc Normanll2007[ ). were responsible for the initial 
metal enrichment of the intergalactic medium (IGM). 
The first supernova (SN) explosions rapidly dispersed 
the heavy elements that were produced during the 
brief lifetime of a Pop III star into the IGM, thus 
beginning the long nucleosynthctic build-up from a 
pure H/He uni verse to one with ubiquit ous metal 
enrichment (e g [Ferrara et al.ll2 000: Madau et al.ll2001l: 
Mori et al."2002'; 'Scannapieco et al. 2002; Thacker et al.' 
2002; Bromm et al. 2003: _Purlanctto fc Locb 2003, 

2004, 



Omukail 120001: iBromm et al. 




-It ^P) \c.f^., 

2001al: IBromm fc Loebl 


20031: ISchneider et all 20031 


20061: iFrebel et all 120071: 



Smith fc Sig urdssoni 120071 ). It is then crucially impor- 
tant to understand the topology of early metal enrich- 
ment and when a certa in region in the un i verse becomes 
supercritical (e. g., Schneider et al.l 120021; iMackey et aU 
20031: IScannapie co ct al. 2003; Ricotti fc Ostrikcr 2004; 
Furlan etto fc Locb 2005; Mattcucc i fc Calura ,2005 



33l: iDaigne et al.l 
~ lYoshida et al.l 



2004 
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Grcif fc Bromm 2006|; |\^nkatcsan 200^). Second, the SN 
blast waves mechanically impacted the halos that hosted 
Pop III stars by heating and subsequently evacuating the 
dense gas inside them. Such a negative feedback effect, 
limiting the capacity for further star formation, has pre- 
viously been considered for low-mass gala xies with cor- 
respondingly shallow po tential w ells fe.g.. 'La rsonlll97^ 
[Dekcl&Silk 1986; Ma c Low fc F errara 19991 but also 
for halos in the vicinity of a S N progenitor, where th e 
net effect is still uncertain (e.g.. lCen fc Riquelmell2007l) . 
A qualitatively different mechanical feedback effect has 
recently been suggested to occur in t he dense, post-shock 
regions of the energet i c blast wave (iMackev et al.l[200l 
iSalvaterra et al.l 120041 : iMachida et all l2005f ) where sec- 
ondary star formation might be induced by the onset 
of gravitational instabilities, thus giving rise to positive 
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feedback. 

The investigation of energetic explosions in the high- 
redshift universe has a long and venerable history. Most 
of this early work has focused on the cooling an d frag- 
mcntation of astrophysical blastwavcs fc.g.. lBertschingeri 
[TqSS: Vishniac et al'l985; W andel 1985,) and th e possi- 
bility of self- propagating gala.xy formation (e. g., Ilkeuchil 
[19811: lOstriker fc CowidflQSl ICarr et al.lll984D . Such ex- 
plosive galaxy formation models typically implied sub- 
stantial distortions to the cosmic microwave background 
blackbody spectrum and were thus e xcluded by the Cos - 
mic Background Explorer {COBE; Fixsen et al.l 119961 ). 
Many of the physical effects discussed in those seminal 
papers, however, remain relevant for our current studies. 

What kind of SN explosion is expected to end the 
life of a massive Pop III star? According to the pre- 
cise progenitor mass, the SN could be either of the con- 
ventional core-collapse type (for masses < 40 Mq) or it 
could be a pair-instability supe rnova (PISN) for masses 
in the range ~ 140-260 Mq ()Heger fc WooslevI [20021: 
iHeger et al.l 120031) . In this paper, we specifically con- 
sider a PISN with a progenitor mass of 200 Mq and 
an explosion energy of Esn = 10^^ ergs. A PISN is 
predicted to completely disrupt the star, so that all 
the heavy elements produced will be released into the 
IGM (e.g.. lFrver et~al][200H[Heger fc Wooslevll2002n . re- 
sulting in extremely high metal yields. However, our 
simulation would also approximately describe a hyper- 
nova explosion, in which a rapi dly rotating , massive star 
undergoes core collap se (e.g., lUmeda fc Nomoto .20021 : 
iTominaga et al.|[2007l ). as far as the energetics and over- 
all dynamics of the blast wave expansion are concerned. 

Theoretically, one would expect that a PISN could 
only have occurred in the early universe, where Pop III 
stars were born massive and mass loss might hav e 
been neghgible (jBaraffe et al.l 120011: IKudritzkil 120021 ). 
The recent discovery of the extremely luminous, rela- 
tively nearby SN 2006gy, which has tentatively been 
interpreted as a PISN (Smith et al. 2007; but see 
Ofek et al. 2007), might defy this theoretical expec- 
tation. Having such a nearby example of a PISN, 
in effect, would provide us with an "existence proof," 
opening up the exciting possibility of detecting the 
first SNe with the upcoming James Webb Space Tele- 
scope (JWST), which will be sensitive e nough to ob- 
serve a single PISN out to z ^ 15 (e.g. , 'Mackey et alj 
2003; Scannapieco et all [20051: IWeinmann & Lilly 200l; 
Wise fc Abe]|l2005t IGardnCT et al. 

In this paper, we carry out cosmological simulations of 
the first SN explosions, using the smoothed particle hy- 
drodynamics code G ADGET in its entro py-conserving 
formulation (ver. 1.1: ISpringel et ani2001[ ). To address 
the questions outlined above, it is crucial to investigate 
the evolution of the SN remnant in three dimensions, 
starting from realistic initial conditions and including 
all the relevant physics. In particular, one needs to im- 
plement a comprehensive model for the cooling and the 
chemical evolution of the g well as an efficient algo- 
rithm to treat the r adiative transfer arou nd the Pop III 
progenitor star (see [Johnson et all 120071 ). Our simula- 
tions greatly improve on earlier work (e.g.. [Bromm et al.l 
[200l in that we calculate the preexplosion situation with 
much greater realism, using a ray-tracing method to de- 
termine the structure and extent of the H ii region. In 



addition, we follow the evolution to much later times, al- 
lowing us to reach the point where the blast wave finally 
stalls and effectively dissolves into the general IGM. We 
are thus able to analyze the expansion and cooling prop- 
erties of the SN remnant in great detail, allowing us to 
draw robust conclusions on the temporal behavior of the 
shock, its morphology, the final shock radius, and the 
total swept-up mass. 

The structure of our work is as follows. In §2, we de- 
scribe the SN progenitor and our method of initializing 
the SN explosion, followed by a test simulation to ver- 
ify the accuracy of our results. Subsequently, we discuss 
the cosmological setup of the main simulation and our 
treatment of the H ii region. We then concentrate on 
the evolution of the SN remnant, finding a simple an- 
alytic model that summarizes its expansion properties 
(§3). In §4, we investigate the mechanical feedback of 
the blast wave on neighboring minihalos and elaborate 
on the prospect of triggering gravitational fragmentation 
in the dense shell. In §5, we discuss the relevance of metal 
cooling, the coarse-grain dispersal of metals, and the gen- 
eral mixing efhciency. Finally, in §6 we summarize our 
results and discuss important cosmological implications. 
For consistency, all distances are physical (proper) unless 
noted otherwise. 

2. NUMERICAL METHODOLOGY 

The treatment of SN explosions in smoothed particle 
hydrodynamics (SPH) simulations is a demanding prob- 
lem, primarily due to the strong discontinuities arising 
at the shock front. One must set up the initial condi- 
tions with great care, and it is necessary to enforce strict 
timescale constraints in calculating the subsequent evo- 
lution. In light of these challenges, we describe the SN 
progenitor and our method of initializing the SN explo- 
sion, followed by a performance test to verify the correct 
behavior of the shock. 

2.1. The Supernova Progenitor 

Although the initial mass function (IMF) of Pop III 
is poorly constrained, numerical simulations have indi- 
cated that stars forming in primordial halos typically 
attain 100 Mq by efficient accretion and might even 
become as massive as 500 (e.g.. lOmukai &: Pallal 
2003; Bromni fc Loebl l200l lO'Shea &: NormanI 120071) 7 
Heger &: WooslevI (j2002i ) have discussed the fate of such 
massive stars and found that in the range 140-260 Mq a 
PISN disrupts the entire progenitor, with explosion en- 
ergies ranging from 10^"'^ to lO^'^ ergs and yields up to 
y~0.5. 

In the present work, we aim to investigate the most 
representative case, assuming a stellar mass of Af, — 
200 Mq. We conservatively adopt an explosion energy 
of Esn — 10^^ ergs and a yield of y — 0.1, to include 
the possibility of a less top-heavy IMF and a hypernova 
explosion, but we caution that such models generally pre- 
dict energy input via bipolar jets, which may invalidate 
the assumption of a spherically symmetric blast wave. 
In light of the uncertain progenitor, we note that the 
dynamics of the SN remnant are governed mainly by the 
explosion energy, together with the IGM density distribu- 
tion, and that the composition and metal content of the 
stellar ejecta become important only at very late times 
when the enriched gas recoUapses (see §5.1). 
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2.2. Energy Injection 

During the early evolution of the SN remnant, i.e., 
typically less than 10'* yr after the explosion, the stellar 
ejecta are confined to a thin shell propagating at constant 
velocity into the surrounding medium, while secondary 
shocks qu ickly bring the interior to a uniform tempera- 
ture (e.g.. IGullll973n . This stage in the evolution of the 
SN remnant is termed the free-expansion (FE) phase, 
and it lasts until the swept-up mass becomes compara- 
ble to the ejecta mass. 

To reproduce these well-known initial conditions, we 
inject the kinetic energy of the SN as thermal energy into 
the iVsn — M^/rusph innermost particles surrounding the 
center of the box, where mgph — 5 Mq corresponds to 
the SPH particle mass. This increases the thermal en- 
ergy of all particles inside rfe to Eth = -Esn/-^sn, and 
the resulting temperature gradient acts as a piston on 
the surrounding gas. The ensuing shock accelerates the 
material in the vicinity of rfc and creates a thin shell of 
highly supersonic material, initiating the Sedov- Taylor 
(ST) phase of the SN remnant. In addition, this method 
allows us to track the dispersal of metals, as each parti- 
cle inside rfc represents the original stellar content. Al- 
though the so-obtained resolution is crude, we can nev- 
ertheless quantify the coarse-grain chemical enrichment 
properties of the SN. 

2.3. Test Simulation 

To test our method of setting up the initial condi- 
tions, and also verify that the code reliably calculates 
the shock propagation, we compare the radial profiles of 
a test simulation with the well-known ST solution. For 
this purpose we place the SN according to the above pre- 
scriptions in a noncosmological box of length 1 kpc with 
200^ gas particles and switch off gravity, chemistry, and 
cooling. The particles are distributed randomly, so that 
the density fluctuates slightly around the mean of nn = 
0.5 cm~^, while the temperature is set to 2 x lO"* K. This 
choice reflects the situation in the vicinity of a 200 Mq 
Pop III star af ter the surrounding medium has been pho- 
toheated (e.g iKitayama et al.ll200 iWhalen et al.l l2004l: 
I Alvarez et alll2006l: I Johnson fc Brommll2007D . but it also 
reproduces the initial conditions of the main simulation 
(see Fig. 8). 

Once we initialize the simulation, the shock propagates 
into the surrounding medium according to the ST solu- 
tion. As the duration and time-dependent properties of 
the gas in this phase are crucial for the late-time behavior 
of the SN remnant, we briefly review the relevant details. 

2.3.1. ST Solution 

For a strong shock in an adiabatic gas, the Rankine- 
Hugoniot jump conditions imply for the density, velocity, 
pressure, and temperature directly behind the shock 
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Fig. 1. — ST solution: density, velocity, pressure, and inverse 
temperature as a function of position relative to the shock. Most 
of the mass piles up just behind the shock and forms a dense shell, 
while the isobaric interior regions remain at high temperatures. 



where p is the density of the surrounding medium, ^ 
its mean molecular weight, mn the mass of the hydro- 
gen atom, and fee Boltzmann's constant. Solving the 
continuity, Euler, and energy equations, one can relate 
these quantities at any point behind the shock to those 
directly behind the shock as a function of position. Fig- 
ure 1 shows the resulting self-similar, time-independent 
character of the ST solution, implying that most of the 
swept-up mass piles up just behind the shock and forms 
a dense shell, while the interior regions remain at high 
temperatures. 

Completing the ST solution, the absolute position of 
the shock as a function of time is given by: 



1/5 



(5) 



where (3 = 1.17 for an adiabatic gas. Thus, given the 
SN explosion energy and the density of the surrounding 
medium, all quantities of relevance around the progenitor 
are well-determined. 

2.3.2. Test Results 

Using the above results, we investigate whether the 
test simulation has successfully reproduced the ST solu- 
tion. Figure 2 compares the profiles of the test simulation 
after 1 Myr, when enough time has passed to enable a re- 
laxation to the asymptotic solution, to the analytic, but 
dimensionalized, ST profiles shown in Figure 1. Apart 
from minor deviations caused by initial density fluctua- 
tions and higher order shocks, in particular in regard to 
the velocity distribution, we find that the simulation and 
the ST solution are in good agreement. The slight offset 
at the position of the shock is inevitable in light of kernel 
smoothing. 

Being assured that we can reliably treat the propaga- 
tion of shocks, we turn to the main focus of this work and 
describe the setup for a SN explosion in the high-redshift 
universe. 

2.4. Main Simulation 
2.4.1. Initial Setup 

To investigate the long-term evolution of the SN rem- 
nant, but also to correctly calculate the structure and 
extent of the preceding build-up of the H ii region 
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Fig. 2. — Density, velocity, pressure, and temperature of the 
shocked gas after 1 Myr. Black dots represent the test simulation, 
while the gray lines show the dimensionalized ST solution. Apart 
from deviations caused by higher order shocks and kernel smooth- 
ing, the simulation reproduces the analytic profiles relatively well. 



around the progenitor, we choose a cosmological box 
of length 150/h kpc (comoving), with 200^^ particles 
per species (DM and gas). We initialize the simu- 
lation at z = 100 deep in the linear regime, adopt- 
ing for this purpose a concordance ACDM cosmology 
with the following parameters: matter density Qm = 
1 — r^A = 0.3, baryon density fib = 0.04, Hubble param- 
eter h ^ Hq/ (lOO km s^^ Mpc^^) = 0.7, spectral index 
1.0, a nd a top-hat fluctuation power erg — 0.9 (e.g. 



iSpergel et a l. 2003). Initial density and velocity pertur- 
bations are imprinted according to a Gaussian random 
field and grow in proportion to the scale factor until the 
onset of nonlinearity. At this point the detailed chemi- 
cal evolution of the gas becomes crucial , and we ap ply 
the same chemical network as in Johnso n et al.l (|2007l ) to 
track the abundances of H, H"*", H~, H2, Hj, He, He"*", 
He"'""'", and e~, as well as the five deuterium species D, 
D+, D~, HD, and HD~. All relevant cooling mechanisms 
in the temperature range 10-10® K are implemented, in- 
cluding H and He resonance processes, bremsstrahlung, 
inverse Compton (IC) scattering, and molecular cooling 
for H2 and HD. Metal cooling does not become impor- 
tant for the entire lifetime of the SN remnant, yet we 
postpone a more detailed discussion of this issue to §5. 
We do not take into account the emission of radiation 
by the post-shock gas, which acts to create a thin layer 
of fully ionized material ahead o f the shock and sup- 
presses molecule forniation (e.s... fShull &: McKed 119791 : 
iShapiro fc Kandll987l : [Kang fc ShaDiralll992D . since (1) 
the SN remnant expands into an H 11 region, and (2) we 
find that molecule formation in the post-shock gas be- 
comes important only at late times, when it has cooled 
to below 10-* K (see §3.4). 

With these ingredients, the first star forms in a halo 
of Mvir ~ 5 X 10^ Mq and ryj^ - 100 pc at z - 20 
in the canonical fashion (e.g.. iBromm "erall[l999l . [200l : 
I Abel et al]|2002f ). We determine its location by identi- 




Redshlft: 
z = 21 

Boxsize: 
150/h kpc [comoving) 

Slice Width: 
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Fig. 3. — Hydrogen number density averaged along the line of 
sight in a slice of 10/ h kpc (comoving) around the first star, form- 
ing in a halo of total mass Mvir ~ 5 X 10^5 Mq at 2 ~ 20. Evidently, 
the host halo is part of a larger conglomeration of less massive mini- 
halos and is subject to the typical bottom-up evolution of structure 
formation. 

fying the first particle that reaches a density of nn = 
10'' cm"'^. At this point the gas "loiters" around a tem- 
perature of 200 K and typically attains a J eans mass of a 
few 10'^ M(T) before further collapsing (e.g. . IBromm et aTl 
[200llGlovej[2005l) . For simplicity, we assume that such 
a clump forms a single star, and we find that its location 
is reasonably well resolved by the minimum resolution 
mass, Mi-OS — 500 Mq. In Figure 3, we show the hydro- 
gen number density in the x-y and y-z plane, centered on 
the formation site of the first star. Evidently, the host 
halo is part of a larger overdensity that will collapse in 
the near future and lead to multiple merger events. This 
behavior is characteristic of bottom-up structure forma- 
tion, and our simulation therefore reflects a cosmological 
environment typical for these redshifts. 

2.4.2. H 11 Region 

The treatment of the H ii region around the star 
is crucial for the early- and late-time behavior of the 
SN remnant. The photoevaporation of the host mini- 
halo greatly reduces the central density and extends the 
energy-conserving ST phase, whereas after an intermedi- 
ate stage the enhanced pressure in the H ii region leads to 
an earlier transition to the final, momentum-conserving 
phase. In addition, the shock fulfills the stalling criterion 
(i.e. fsh = Cs, where Cs is the sound speed of the photo- 
heated IGM) much earlier in the H II region compared to 
previously unheated gas. We have found that neglecting 
the presence of the H ii region around the star, which 
extends well into the IGM, leads to a final shock radius 
a factor of 2 larger, demonstrating its importance for the 
long-term evolution of the SN remnant. 

To determine the size an d structure of the H i i region, 
we proceed analogously to iJohnson et al.l ()2007[ ). In de- 
tail, we initially photoheat and photoionize a spherically 
symmetric region surrounding the star up to a maximum 
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distance of 200 pc, where we find a neighboring mini- 
halo. We determine the necessary heating and ioniza- 
tion rates by usin g the properties of a 200 Mp^ Pop II I 
star found by Bro mm et al.l (l2001bn and lSchaere^ (|2002f ). 
After about 2 Myr, when the star has reached the end 
of its lifetime, the hydrodynamic shock has propagated 
to rvir/2 and photoevaporated the central regions of the 
host halo. Figure 8a shows the resulting density, tem- 
perature, pressure, and velocity profiles, which d isplay 
the characteristics of the analvtic lShu et ahl (j2002f ) solu- 
tion (i.e. pressure equilibrium throughout the interior, 
while the density [temperature] become almost constant 
at small [large] radii). In our case, the average inte- 
rior density drops to nu — 0.5 cm^'^, while the cen- 
tral temperatures rise to roughly 4 x lO** K (see Fig. 
8 a). The radial profiles agree relatively well with the 
self-consis tent radiation-hyd rodyn amics simulations per - 
formed in lAbel et all (|2007D and lYoshida et all (|2007D . 
with slight differences most likely due to the harder spec- 
trum and shorter lifeti me of a 200 M(7) star comp ared to 
a 100 A/q star (see also lJohnson fc Bromml[2007[ ). 

Since ionizing radiation escap es the minihalo a fter only 
a few thousand years (e.g., Al varez et al.|[2006l ). we de- 
termine the final size and structure of the H ii region 
by performing the r ay-tracing algorithm introduced in 
iJohnson et al.l (|2007| ). This routine finds the Stromgren 
radius along 10^ rays around the star, with each ray con- 
sisting of 200 radial bins. For this purpose we assume 
that helium has the same ionization properties as hydro- 
gen, and we do not treat a separate He ill front. Once an 
individual cell fulfills the Stromgren criterion, we ionize 
its content and set its temperature and molecule abun- 
dances to the values at the outer edge of the photoheated 
region. This method is appropriate for the propagation 
of ionization fronts in the general IGM, but it does not 
correctly treat the photoevaporation of neighboring mini- 
halos, when the ioni zation front becomes D type (e.g., 
lAhn k ShaT)irol[200l . This issue is particularly impor- 
tant with respect to feedback on neighboring minihalos 
and is further discussed in §4. 

Figure 4 shows the gas temperature after the main- 
sequence lifetime of the star, and thus indirectly the size 
and structure of th e H ii region. In agreement with 
lAlvarez et al.l (j2006D . we find that the H ii region can 
be as large as 5 kpc, but efficient shielding by neighbor- 
ing minihalos may limit its size to a few 100 pc in some 
directions. The H ii region initially cools via IC scat- 
tering, while the pressure gradient at the boundary of 
the H II region leads to a gradual adiabatic expansion on 
timcscales comparable to the Hubble time at z ~ 20 (see 
Figs. 4 and 9; also see Johnson & Bromm 2007). Even 
though molecule fractions rise to typically x-^^ ~ 10~^ 
and xhd ~ 10~^, molecular cooling remains inefficient 
due to the low densities of the photoheated gas (e.g., 
IJohnson k, Bro mm 2007) . We do not treat the evolution 
of a separate Lyman- Werner (LW) front, as molecules 
inside the H ii region are destroyed primarily by coUi- 
sional dissociation and charge transfer, while they are 
quickly reformed in more massive, neigh boring miniha- 
los o i ice the central sourc e turns off (see IJohnson et all 
[20071 lYoshida et al.l[200l . 




Fig. 4. — Temperature averaged along the line of sight in a slice of 
10 /h kpc (comoving) around the star after its main-sequence life- 
time of 2 Myr. Ionizing radiation has penetrated nearby minihalos 
and extends up to 5 kpc around the source, heating the IGM to 
roughly 2 X 10 K, while some high-density regions have effectively 
shielded themselves. 

In the course of the simulation, neighboring miniha- 
los might reach high enough densities to form stars, and 
the subsequent photoheating could significantly infiuence 
the propagation of the SN remnant. In the present work, 
we aim to investigate a single unperturbed SN explosion, 
and thus we rule out further star formation by employing 
the sink particle algorithm used in Johnson et al. (2007). 
This routine forms sink particles once the hydrogen num- 
ber density exceeds 10"* cm^^ and further accretes par- 
ticles inside the Bondi radius (|Bondilll952D . This proce- 
dure prevents disturbances arising from the expansion of 
additional ionization fronts and allows us to concentrate 
on the feedback caused by the SN explosion. 

With these preparations in place, we reinitialize the 
simulation at the end of the EE phase according to §2.2. 
In the following we discuss the evolution of the SN rem- 
nant until it effectively dissolves into the IGM. 

3. EXPANSION AND COOLING PROPERTIES 

Throughout its lifetime, the remnant goes through four 
evolutionary stages, each of which is characterized by a 
different physical mechanism becoming dominant. Based 
on this chronological sequence, we discuss the expansion 
and cooling properties of the SN remnant with respect to 
the simulation results, and we summarize the remnant's 
behavior with a simple analytic model. 

3.1. Phase I: Free Expansion 

At very early times, the SN remnant enters the EE 
phase and propagates nearly unhindered into the sur- 
rounding medium. It expands with a constant velocity 
of v^^ = 2Esn/Mej, where for our case M^j — A'E. The 
duration of the EE phase is given by tfc = rtc/uoj, or 



2.4.3. Sink Particles 



(6) 
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where rte is the radius at which the swept-up mass equals 
the mass of the original ejecta; i.e, 



\47rmH'T-H 



1/3 



(7) 



where X = 0.76 is the pr imordial mass fraction of hydro- 
gen (e.g.. lKitavama fc Yoshida 20051. 

After tfe, the inertia of the swept-up mass becomes 
important, and the shock undergoes a transition to the 
ST phase. Since we do not explicitly model the FE phase, 
we use the above analytic arguments to obtain rte < 
20 pc and < lO"* yr (see Fig. 7). 

3.2. Phase II: Sedov- Taylor Blast Wave 

According to §2.2, the simulation begins when the 
shock undergoes a transition to the ST phase. Since 
we have already discussed the properties of the ST solu- 
tion in §2.3, we can apply these results to understand the 
early behavior of the SN remnant. With the knowledge 
that the previous photoheating has created an average 
density profile of nn — 0.5 cm~^ in the vicinity of the 
progenitor (Fig. 8a), we find with the help of equation 
(5) that the shock approaches rvir/2 after about 10^ yr. 
At this point it catches up with the previously estab- 
lished photoheating shock, where the outlying density 
profile becomes isothermal. The simulation results (Fig. 
86) agree with this analytic prediction, yet the profiles 
do not quite resemble the ST solution, partly because 
not enough time has passed, but also because cooling 
in the dense shell becomes important, violating energy 
conservation. Consequently, the ST phase ends and the 
remnant undergoes a second transition. 

To understand the origin of this transition, we elabo- 
rate on the cooling mechanisms responsible for radiating 
away the thermal energy of the remnant. At the high 
temperatures behind the shock, these are H and He res- 
onance processes (RPs), bremsstrahlung (BS), and IC 
scattering, with the former consisting of coUisional ion- 
ization, excitation, and recombination cooling. Figure 
5 shows the relevant cooling rates, normalized to a hy- 
drogen number density of unity, where the abundances 
of the various species follow from assuming coUisional 
ionization equilibrium. Evidently, IC scattering and BS 
are important at high temperatures, whereas RPs be- 
come significant below 10^ K. However, BS and RPs are 
proportional to the density squared, while IC scattering 
exhibits only a linear dependence. Thus, RPs will be 
the most important coolant in the dense shell, while IC 
scattering is dominant in the hot interior. 

At z = 20, the cooling time for IC scattering is ap- 
proximately 10 Myr, independent of temperature and 
density, while cooling times for RPs in the dense shell 
are generally much lower. To determine when radia- 
tive losses affect the energetics of the SN remnant, we 
must equate the shell cooling time with the expansion 
time. Figure 6 shows this analytic relation for the initial 
conditions of the simulation, i.e. Esn — 10^^ ergs and 
TT-H — 0.5 cm""^. Confirming the simulation results, this 
prediction yields that RPs efficiently cool the dense shell 
to lO"* K after about 10^ yr. At this point the shocked 
gas separates into a hot, interior bubble with tempera- 
tures above 10^ K and a dense shell at lO** K bounded 
by a high-pressure gradient. This multi-phase structure 




Fig. 5. — Cooling rates as a function of temperature at 2 ~ 20 
for H and He RPs (solid line), BS (dotted line), and IC scattering 
(dashed line) in coUisional ionization equilibrium, with njj set to 
unity. At low temperatures RPs dominate, while above 10^ K IC 
scattering and BS become important. However, RP's and BS are 
proportional to the density squared, while IC scattering exhibits 
only a linear dependence. 




Fig. 6. — Cooling time over expansion time in the dense shell 
for RPs (solid line), BS (dotted line), and IC (dashed line). The 
long-dashed line indicates when cooling becomes efficient, implying 
that RPs are responsible for ending the ST phase after 10"" yr. 
At z = 20, the cooling time for IC scattering is roughly 10 Myr, 
independent of temperature and density. 



is clearly visible in Figure 8&, which remains intact for 
< 10 Myr, when IC scattering becomes important and 
cools the last remnants of the interior bubble to the tem- 
perature of the dense shell. With energy conservation 
no longer valid, the ST phase ends before the shock has 
relaxed to the asymptotic solution, yet we nevertheless 

2/5 

find that Vgh oc t^^ fits the temporal scaling of the mass- 
weighted mean shock radius relatively well (see Fig. 7). 



3.3. Phase III: Pressure- Driven Snowplow 

After the ST phase ends, the pressurized, interior bub- 
ble drives a dense shell, and one speaks of a pressure- 
driven snowplow (PDS). To analytically describe the fur- 
ther evolution of the SN remnant, we assume that the en- 
tire swept-up mass Mgw is confined to an infinitely thin 
shell. In light of the steep density profile toward the inte- 
rior, this assumption is justified, leading to an equation 
of motion of the form: 



d {MswVs 
dt 



AwrlPb 



(8) 
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where Pb is the pressure of the hot, interior bub- 
ble, and the external pr essure has been neglected (e.g., 
lOstriker fc McKed [Tosl . Since Msw esc rsh in an r 
density profile and oc in the adiabatically ex- 
panding interior, one can solve the above equation with 

2/5 

a power law of the form r^h oc t^^ . Interestingly, this 
procedure yields the same temporal scaling as the ST so- 
lution, and we therefore do not expect a change in slopes 
after 10^ yr. This is a direct consequence of the transition 
to the PDS phase once the shock approaches the isother- 
mal density profile in the outskirts of the halo. Figure 7 
confirms this result, as the simulation shows only a slight 

2/5 

deviation from the analytically derived t^^ slope. 

Using the above model, we can further determine when 
the PDS phase ends. The pressure directly behind the 
shock after 10^ yr can be estimated with equation (3) as 
Pb/fce — 3 X 10^ K cm~^, in agreement with the simula- 
tion (see Fig. 86). With Pb oc and Vsh oa t^^ , we fur- 
ther find Pb oc , implying that after roughly 1 Myr the 
interior pressure has dropped to Pb/fce — 3 x 10"* K cm~^ . 
At this point pressure equilibrium between the hot inte- 
rior and the dense shell has been established, and the 
shock is driven solely by its accumulated inertia. Fig- 
ure 8c confirms this prediction, showing that the interior 
pressure has indeed dropped to that of the dense shell. 
In contrast, the temperature only drops proportional to 
^sh and remains high. 

3.4. Phase IV: Momentum-Conserving Snowplow 

With the pressure gradient no longer dominant, the 
SN remnant is driven by the accumulated inertia of the 
dense shell and becomes a momentum-conserving snow- 
plow (MCS). In analogy to the derivation performed in 
§3.3, the position of the shock as a function of time can 
be obtained by solving equation (8) in the absence of a 
pressure term. Since the shock has not yet propagated 
beyond the surrounding density profile, this yields 

an initial scaling of rgh oc t^'^ . 

At later times, the shock finally leaves the host halo 
and encounters neighboring minihalos in the y-z plane, 
but underdense voids perpendicular to the y-z plane (see 
Fig. 3). This increases its radial dispersion, while at the 
same time the Hubble expansion becomes important and 
serves to expand the medium on which the SN remnant 
propagates, thus raising its physical shock velocity. We 
can therefore only estimate the temporal scaling of the 
mass-weighted mean shock radius based on the simula- 
tion results, finding, interestingly, that rgh maintains its 

1/2 

scaling until it fulfills the stalling criterion (see Fig. 
7). This occurs after about 200 Myr, close to the Hubble 
time at z ~ 20, when the shock velocity approaches the 
local IGM sound speed and becomes indistinguishable 
from sound waves. Evidently, the increased slope at late 
times is a fundamental difference between SNe occurring 
at high redshifts and in the present-day universe, where 

^sh OC tj^ . ^ 

In the meantime, the hot interior has expanded adia- 
batically for roughly 10 Myr, when IC scattering becomes 
important, quickly cooling the interior to temperatures 
of the dense shell (see Fig. 8d). Even later, the high elec- 
tron fraction persisting in the interior and the dense shell 



leads to efficient molecule formation, and typical abun- 
dances of ^ 10~^ and xhd ^ 10^'' are reformed. 
By the end of the simulation, both phases have adiabat- 
ically cooled to T ~ 10'^ K, while densities in the dense 
shell have dropped to nn ~ 10^^'^ cm~'^ (see Figs. 9 
and 14). The interior is even more underdense, and we 
find that the SN explosion has completely disrupted the 
host halo, preventing further star formation inside of it 
for at least a Hubble time at z ~ 20. Due to the continu- 
ous adiabatic expansion of the post-shock gas, molecular 
cooling remains inefficient for the entire lifetime of the 
SN remnant. 

3.5. Summary of Expansion Properties 

Summarizing the expansion properties of the SN rem- 
nant. Figure 7 shows the mass- weighted mean shock ra- 
dius as a function of time, together with the analytically 
derived power laws. The good agreement indicates that 
one can quantify the evolution of the SN remnant by 
means of simple physical arguments and obtain relatively 
accurate results. 

We find a final mass-weighted mean shock radius of 
2.5 kpc, which is roughly a factor of 2 smaller than the 
H II region. Figure 7 also shows the radial dispersion of 
the shock, indicating that the surrounding medium be- 
comes highly anisotropic once the SN remnant leaves the 
host halo. After about 5 Myr, the shock encounters the 
first neighboring minihalos in the y-z plane and begins 
to stall, while it continues to propagate into the voids 
perpendicular to the y-z plane. The vast majority of 
the shocked material resides in the general IGM, as the 
mass-weighted mean closely traces the maximum shock 
radius. 

To conclude this section. Figure 10 shows the swept- 
up gas mass as a function of time for the simulation and 
the analytic model, where the slopes have been deter- 
mined by the relation Msw oc pr^^, except at very late 
times when the shock structure becomes too complex 
for a simple analysis. We find a final swept-up mass of 
2.5 X 10^ Mq, which is roughly a factor of 2 smaller than 
the mass enclosed inside a sphere of radius 2.5 kpc, as- 
suming the background density at z ~ 12. This further 
demonstrates the overall asymmetry of the shock (see 
Fig. 14). 

4. FEEDBACK ON NEIGHBORING HALOS 

So far we have concentrated on the propagation of the 
SN remnant into the IGM and treated neighboring mini- 
halos as disturbances to the radial density profile. In the 
following, we explicitly investigate the mechanical feed- 
back of the SN on nearby minihalos and discuss possible 
consequences for star formation. For this purpose we 
have carried out two additional simulations, one without 
feedback and another with radiative feedback, which en- 
able us to disentangle the effects of photoheating and the 
SN shock. 

Combining the so-obtained results. Figure 11 shows the 
distances of all star-forming halos from the initial SN pro- 
genitor as a function of expansion time, when they have 
reached the threshold density nn = 10"* cm~'^. Following 
our argumentation in §2.4, the halo irreversibly collapses 
at this point and forms a star. The shades of the sym- 
bols in Figure 11 indicate their affiliation to the different 
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Fig. 7. — Evolution of the SN remnant in its entirety, from z ~ 20, 
when the SN explodes, to 2 ~ 12, when the shock finally stalls. Its 
total lifetime is about 200 Myr, or a Hubble time at z ~ 20. The 
black dots indicate the mass-weighted mean shock radius accord- 
ing to the simulation, while the dashed line shows the analytic 
solution. For both we find a final mass-weighted mean shock ra- 
dius of 2.5 kpc. The shaded region shows the radial dispersion 
of the shock, indicating that it increases significantly once the SN 
remnant leaves the host halo and encounters the first neighboring 
minihalos in the y-z plane. The bulk of the SN remnant propa- 
gates into the IGM, since the mass- weighted mean closely traces 
the maximum shock radius. 

simulations (i.e. black, dark gray, and light gray sym- 
bols represent the no-feedback, photoheating-only, and 
main simulation runs), while the shapes of the symbols 
represent individual halos. For orientation, we show the 
mass- weighted mean shock radius at late times according 
to Figure 7. 

4.1. Delay by Photoheating 

In the photoheating-only case, we find that the col- 
lapse of all three star-forming halos is delayed (see Fig. 
11). While we do not properly resolve D-type ioniza- 
tion fronts that may develop within shielded minihalos, 
we find that ionizing radiation can penetrate deep into 
their cores and suppress cooling and the accretion of gas. 
The distance, mass, and maximum density of the halo at 
the onset of photoheating are all crucial for the extent 
of the delay. The most massive halo in our simulation, 
halo a, experiences a delay of only 25 Myr, while halos 
b and c are less dense at the time of photoheating and 
experience delays in excess of 80 Myr. This is rou ghly 
consistent with the results of lAhn fc Shapiro! (|2007D who 
argue that radiative feedback on neighboring minihalos 
strongly depen ds on their evolutionary stag e (see also 
iMesinger et 311120061 : Isiisa fc Umemurall2006f i. For more 
quantitative results, one must perform self-consistent 
radiation-hydrodynamics simulations that include the ef- 
fects of photodissociating radiation fe.g.. lYoshida et all 
[200l . 

Due to the inaccuracies mentioned above, we cannot 
draw robust conclusions on the subsequent mode of pri- 
mordial star formation, i.e., Pop III or Pop II. 5, but 



we anticipate this to be a function of th e state of the 
collapse and the ionizing fl u x fsee also iMackey et al.l 
200f 'Greif fc Bromml [200l: iJohnson fc Bromml l200a 
Jappsen et al. 2007.) . 

4.2. Shock-driven Collapse 

In contrast to feedback by ionizing radiation, we find 
that in our case the shock of the SN remnant acts to en- 
hance halo collapse and promote star formation. This be- 
havior is evident from Figure 1 1 , as halos a and b collapse 
about 15 Myr earlier with respect to the photoheating- 
only run. In both cases, the shock compresses the dense 
cores of each of the halos, which serves to enhance cool- 
ing and accretion. Although our resolution is too crude 
for quantitative conclusions, wc find that the mechanical 
impact of the SN remnant mitigates the delay caused by 
photoheating and leads to a slightly increased star for- 
mation rate. This situation might be different for star- 
forming halos much closer to the SN explosion, where 
ram pressure stripping could be sufficient to dispel a large 
fraction of the gas. However, in the present simulation 
we find that most halos massive enough to form stars are 
at sufficient distances from the SN progenitor. 

We note that we neglect the possible effects of pho- 
todissociating radiation emanating from the relic H ii 
region, as it has been shown that recombination radia- 
tion can only produce a relatively weak photodissociating 
flux, which quickly dies away as the gas in the relic H ii 
region recombines on timescales of the order of 1 Myr 
(e.g., Johnson & Bromm 2007; Yoshida et al. 2007; 
see also Glover 2007). In addition, we are able to ne- 
glect the radiation generated by the supernova-shocked 
gas as a source of photodissociating radiation concern- 
ing neighboring minihalos, since the shock is only mov- 
ing at velocities < 20 km s~^ when it arrives at the 
star-forming minihalos mentioned above. This implies 
a photodissociating flux of Jlw ^ 10^^, where Jlw is 
in un its of IQ-^i ergs s"! cm^^ Hz"^ sr"! ()Shull fc Silkl 
Il979f l. Such a low level of flux will not have a substantial 
impact on the evolution of gas in collapsing minihalos, 
as it implies a timescale for photodissociation of H2 of 
the order of 1 Gyr, much longer t han the Hubble time 
at the rcdshifts wc consider (e.g., lOh fc HaimanI [20031 : 
Ijohnson fc Bromm 2006). However, a minihalo within 
~ 200 pc of the SN explosion may be subject to a signif- 
icant photodissociating flux, owing to the higher shock 
velocity at small distances. 

4.3. Mixing Efficiency 

An important question is whether a large fraction of 
metals can penetrate neighboring halos and efficiently 
mix with their cores, thus changing the mode of star 
formation from Pop Ill/Pop II. 5 to Pop II. Although 
a detailed quantitative analysis would require dedicated 
high-resolution simulations that are unavailable here, we 
can nevertheless estimate the mixing efficiency by apply- 
ing the criterion for the operation of Kelvin-Helmholtz 
(KH) instabilities given by Murray et al. (1993; see also 
Wyithe fc Cen 2007; Gen fc Riquelme 2007): 
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Fig. 8. — Density, velocity, pressure, and temperature around the SN progenitor star immediately before its death, and 10^, lO", and 
10^ yr after the SN explosion. Black dots represent normal SPH particles, while gray dots represent the initial stellar ejecta. (a): The 
hydrodynamic shock of the H II region has approache d r^„l2 at ^_50 pc after establishing pressure equilibrium in the interior, while the 
density and velocity profiles assume the characteristic IShu et all 1)20021 ) solution of a champagne flow. The temperature of the H II region 
is typically 2 X 10* K. (h): The shock created by the SN explosion propagates into the surrounding medium according to the ST solution, 
with the characteristic formation of a dense shell and a hot, interior bubble. After 10^ yr, the shock approaches the previous photoheating 
shock with expansion velocities well in excess of 100 km s~^. At the same time, cooling in the dense shell by RPs becomes efficient and 
the ST phase ends. The resulting sharp temperature drop between the interior bubble at T > 10^ K and the dense shell at T ~ 10* K is 
clearly visible, (c): The SN remnant has undergone a transition to the PDS phase and is driven by the high pressure in the interior and the 
momentum of the dense shell. This phase ends after about 1 Myr, when pressure equilibrium has been established. The pressure drop in 
comparison with (h) is clearly visible, while temperatures in the interior remain high, (d): After 10 Myr, the radial dispersion of the shock 
has increased dramatically, ranging from > 200 pc where the shock encounters the first neighboring minihalos, to < 1 kpc perpendicular 
to the overdensities in the y-z plane. 
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Fig. 9. — Temperature averaged along the line of sight in a slice of 10/h kpc (comoving) around the x-y plane after 1, 10, 50, and 200 Myr. 
In all four panels, the H II region and SN shock are clearly distinguishable, with the former occupying almost the entire simulation box, 
while the latter is confined to the central regions, (a): The SN remnant has just left the host halo, but temperatures in the interior are 
still well above 10* K. (b): After 10 Myr, the asymmetry of the SN shock becomes visible, while most of the interior has cooled to well 
below 10* K. (c): The further evolution of the shocked gas is governed by adiabatic expansion, and its morphology develops a "finger-like" 
structure (see also Fig. 13). (d): After 200 Myr, the shock velocity approaches the local sound speed and the SN remnant stalls. By this 
time the post-shock regions have cooled to roughly 10^ K. 
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Fig. 10. — Swept-up gas mass as a function of time for the sim- 
ulation (black dots) and the analytic model (dashed lines), where 
the latter is related to the shock radius by Afaw oc pr^^^, except 
at very late times, when the shock structure becomes too complex 
for a simple analysis. For both we find a final swept-up mass of 
2.5 X 10^ Mq. 

where g is the gravitational acceleration at the virial ra- 
dius, and D the density ratio of gas in the halo compared 
to the dense shell. For the nearest star-forming halo with 
about 5 X 10^ Mq and rvir — 100 pc, we find I? ~ 10 and 
rsh — 20 km s~^. The left-hand side of equation (9) is 
thus of order 1/10, suggesting that the gas at the outer 
edge of the halo will be disrupted and mix with the en- 
riched material in the dense shell, while the core of the 
halo will remain pristine and stable. This is consistent 
with th e results o f the d etailed simulations conducted by 
iCen fc Riauelmi (|2007l ). and by the same arguments we 
find that halos b and c remain pristine as well. 

Once again, we note that massive minihalos very close 
to the SN progenitor might experience a different behav- 
ior, as expansion velocities are of the order 100 km 
when the SN remnant leaves the host halo (see Fig. 8). 
In such cases, the shock could disrupt the dense cores and 
trigger efhcient mixing, leading to Pop II star formation 
once the gas recollapses. 

4.4. Gravitational Fragmentation 

Various authors have suggested that cooling in the 
dense shell might lead to gravitational fragmentation 
and trigger secondary star formation (e.g.. | M ackev et al.l 
[200l iSalvaterra et al.l[200l iMachida et alJ l2005f. How- 
ever, due to the previous photoheating, the density of 
the surrounding medium is sufficiently lowered such that 
molecule formation is initially inefficient. Even though 
we do not fully resolve the radiative shock front, we thus 
find that the continuous adiabatic expansion of the dense 
shell renders molecular cooling unimportant and does not 
trigger gravitational instabilities (sec §3.4 and Fig. 8). 
We note that the potential mixing of metals into the 
dense shell and the associated additional cooling might 
affect these results but postpone a more detailed discus- 
sion of this issue to future work. 

Another scenario for secondary star formation is at the 
interface of colliding SN remnants. If such encounters 
take place in significantly overdense regions and occur 
relatively early in the evolution of individual SN rem- 
nants, densities might become high enough for gravita- 
tional fragmentation to occur. 

5. CHEMICAL ENRICHMENT 
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Fig. 11. — Collapse times and distances from the SN progenitor 
for all star-forming minihalos affected by the SN shock. The shades 
of the symbols indicate their affiliation, i.e. black, dark gray, and 
light gray symbols represent the no-feedback, photoheating-only, 
and main simulation runs, respectively, while the shapes of the 
symbols denote the individual halos. For orientation, the dashed 
line shows the mass-weighted mean shock radius at late times ac- 
cording to Fig. 7. In our case, photoheating significantly delays 
star formation, while the SN shock acts to compress gas in neigh- 
boring minihalos and slightly accelerates their collapse. 

We have previously argued that mixing of enriched ma- 
terial with gas in existing star-forming halos is generally 
inefhcient (see §4.3), indicating that the dispersal of met- 
als can only occur via expulsion into the IGM. In §3.5 we 
have shown that the bulk of the shock propagates into 
the voids surrounding the host halo, and we expect that 
chemical enrichm ent proceeds via the same channel (e.g., 
IPieri et aTl 120071 ). This realization is crucial, as the de- 
tailed distribution of metals not only governs the transi- 
tion to secondary star formation, but also the properties 
of the first galaxies, when the enriched material recol- 
lapses into larger potential wells. 

We shed some light on these issues by determining the 
influence of metal cooling on the evolution of the SN 
remnant and by discussing the transport of metals into 
the IGM. Although our resolution is too crude for a de- 
tailed analysis, we can nevertheless discuss the mixing 
efficiency in a qualitative manner and draw some pre- 
liminary conclusions on the ultimate fate of the expelled 
metals. 

5.1. Heat Conduction and Metal Cooling 

During the first 10 Myr, the initial stellar ejecta ex- 
pand adiabatically and apparently do not mix with the 
surrounding material (see Fig. 8). In reality, however, 
the high electron mean free path behind the shock leads 
to heat conduction and gas from the d ense shell evapo- 
rates into the hot, interior bubble fe.g.. lGulil ll973l. This 
effect is not included in the simulation, as it would only 
marginally affect the dynamics of the SN remnant, but 
the onset of Rayleigh- Taylor (RT) and KH instabilities 
is predicted to efficiently mix the metals with primor- 
dial material evaporated from the dense shell, reducing 
the metallicity of the interior by a factor of a few to at 
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Fig. 12. — Cooling time vs. expansion time in the interior. The 
dotted hne represents cooUng by metal lines for a total metallicity 
of Z = Zq, but consisting solely of C, O, Fe, and Si (in equal 
parts), while the dashed line represents IC cooling. The long- 
dashed line indicates when cooling becomes efficient, implying that 
metal cooling becomes momentarily important after 10 Myr, when 
the interior temperatures have dropped to 10^ K and cooling rates 
peak. However, IC scattering becomes efficient at the same time 
and rapidly cools the enriched regions to below 10* K, where metal 
cooling rates drop by a few orders of magnitude. Consequently, the 
presence of metals does not effect the dynamical evolution of the 
SN remnant. 

most one order of magnitude (e.g.. iMadau et aT]|2001[ ). 
In the following, we assume that the overall metallicity 
has dropped by a factor of 5, to Z — Zq, in the course 
of the first million years. 

To determine the importance of metal cooling, we pro- 
ceed analogously to §3.2 and compare the expansion and 
cooling timescales of the interior bubble. F or this pur- 
pose w e adopt the cooling rates provided in iMaio et all 
(|200l for gas in coUisional ionization equilibrium en- 
riched to Z = ^0 but consisting solely of C, O, Fe, and 
Si (in equal parts). Although this does not represent the 
specific yield of a PISN or a hypernova, we here only in- 
tend to give an approximate argument and defer a more 
precise treatment to future work. 

Assuming the time-dependent properties of the sim- 
ulation, we find that metal cooling is only briefly im- 
portant after 10 Myr (see Fig. 12), when the interior 
has adiabatically cooled to 10^ K and metal cooling is 
most efficient. At the same time, however, IC scatter- 
ing becomes important and temperatures quickly fall to 
10^ K, where cooling rates drop by a few orders of mag- 
nitude. Due to the low densities in the interior, the onset 
of fine structure cooling below 10^ K also proves ineffi- 
cient, and we conclude that even for initial metal yields 
as high as ?/ = 0.1, metal cooling is unimportant for the 
entire dynamical evolution of the SN remnant. We em- 
phasize that the presence of metals will become crucially 
important once the enriched gas has recoUapsed into a 
sufficiently massive halo later on, thereby reaching high 
densities again. 

5.2. Instabilities and Distribution of Metals 

Due to inefficient cooling, the evolution of the metal- 
enriched, interior bubble is governed by adiabatic expan- 
sion and preferentially propagates into the cavities cre- 
ated by the shock (see Figs. 8 and 14). Once the shock 
leaves the host halo and becomes highly anisotropic, the 
interior adopts the same behavior and expands into the 
voids surrounding the y-z plane in the shape of an hour- 



glass, with a maximum extent similar to the final mass- 
weighted mean shock radius. This behavior is evident in 
Figure 14, where we plot the hydrogen number density 
and distribution of metal particles at various times after 
the SN explosion. Furthermore, since the radius of each 
metal particle is determined by the smoothing length, we 
find that at late times the interior becomes substantially 
mixed with the initial stellar ejecta. 

When the shock finally stalls, the interior bubble is in 
pressure equilibrium with its surroundings, but it stays 
confined within the dense shell. To investigate the im- 
portance of RT instabilities in this configuration, we esti- 
mate the mixing length Art for large density contrasts be - 
tween two media according to fe.g.. IMadau et alll2001f) : 

Art =^ 27rg4 , (10) 

where g is the gravitational acceleration of the host halo. 
Estimating the distance of the dense shell from the host 
halo with the final mass-weighted mean shock radius, we 
find that mixing between the dense shell and the inte- 
rior bubble takes place on scales < 10 pc in the course 
of a few times 10 Myr. In light of the substantial extent 
of the interior, we conclude that such mixing is gener- 
ally inefficient and that much larger potential wells must 
be assembled to recollect and mix all components of the 
shocked gas. Specifically, turbulence arising in the viri- 
alization of t he first gala xies could be an agent for this 
process fe.g.- IWise fcAbel 2007 ). 

On the other hand, mixing could be important in the 
y-z plane, where, due to the progression of structure 
formation, material falls in along filaments toward the 
metal-enriched gas near the host halo. Since this mech- 
anism will only affect a limited volume, while most met- 
als are expelled into the voids (see Fig. 14), we again 
conclude that a fundamental transition in star formation 
from Pop Ill/Pop II. 5 to Pop II requires the assembly of 
a DM halo with a sufficiently deep potential well. 

To find the minimum mass necessary to recollect the 
hot and underdense post-shock gas residing at T ~ 10"^ K 
and nn ~ 10~^'^ cm~'^ (see §3.4), we estimate by means 
of the cosmological Jeans criterion that a DM halo of at 
least Mvir ^ 10* Mq must be assembled. Assuming that 
in the course of its virialization the initial stellar ejecta 
mix with the entire swept-up mass, we find a final, av- 
erage metallicity of Z ~ 10^^-^ Zq. Although this value 
is well above Zcnt, we emphasize that the final topol- 
ogy of metal enrichment could be highly inhomogeneous, 
with pockets of highly enriched material on the one hand, 
but regions with a largely primordial composition on the 
other hand. 

6. SUMMARY AND CONCLUSIONS 

We have investigated the explosion of a 200 Mq 
PISN in the high-redshift universe by means of three- 
dimensional, cosmological simulations, taking into ac- 
count all necessary chemistry and cooling. Using a ray- 
tracing algorithm to determine the size and structure of 
the H II region around the progenitor star, we have fol- 
lowed the evolution of the SN remnant until it effectively 
dissolves into the IGM, and discussed its expansion and 
cooling properties in great physical detail. Specifically, 
we have found that a chronological sequence in its evo- 
lution, based on various physical mechanisms becoming 
dominant, allows the introduction of a simple analytic 
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Fig. 13. — Three-dimensional view of the SN remnant at z = 15, or 100 Myr after the SN explosion, when it has reached a radius of 
2 kpc. The finger-like morphology of the shock becomes visible as the SN remnant propagates at varying speeds in different directions, 
caused by an anisotropic density distribution. 



model summarizing its expansion properties. The SN 
remnant propagates for a Hubble time at ^ ~ 20 to a fi- 
nal mass- weighted mean shock radius of 2.5 kpc, roughly 
half the size of the H ii region, and in this process sweeps 
up a total gas mass of 2.5 x 10^ Mq. We have found that 
its morphology becomes highly anisotropic due to en- 
counters with filaments and neighboring minihalos in the 
y-z plane, but underdense voids perpendicular to the y-z 
plane. Based on the high explosion energy, the host halo 
is entirely evacuated, while in our case shock compres- 
sion of neighboring minihalos partially offsets the delay in 
star formation due to negative feedback from photoion- 
ization heating. In contrast, we do not observe gravita- 
tional fragmentation triggered by efficient cooling behind 
the SN shock, which could in principle lead to secondary 
star formation. We have found that the metal-enriched, 
interior bubble expands adiabatically into the cavities 
created by the shock and preferentially propagates into 
the voids of the IGM with a maximum extent similar to 
the final mass- weighted mean shock radius. We have es- 
timated that RT instabilities do not efficiently mix the 
dense shell with the interior, but that material falling in 
along filaments could mix with metal-enriched gas near 
the host halo. Finally, we have concluded that a DM 



halo of at least Mvh- ^10® Mq must be assembled to 
recollect and mix all components of the shocked gas. 

Based on the simulation, we find that a single PISN 
can enrich the local IGM to a substantial degree. If en- 
ergetic SNe were indeed a common fate for the first stars, 
they might have deposited metals on large scales before 
massive galaxies formed and outflows were suppressed by 
their increasingly deep potential wells. Hints on ubiqui- 
tous metal enrichment have recently been found in the 
low column density Lyg forest ( e.g., 'So ngaila fc Cowi"i 
1996; Songaila 2001} lAguirre et al. 2005) a nd in dwarf 
spher oidal satellites of the Milky Way (jHelmi et al.l 
I2006f ). where the presence of a "bedrock" metallicity was 
inferred. Based on these observations, various authors 
have argued that SNe occurring in the shallow poten- 
tial wells of mi nihalos could subst a ntiall y preenrich the 
universe (e.g., 'Daiene et al.' '2004', '2006'; 'Yoshida et alj 
[2004; Mattcucci & Calura 2005; Grcif & Bromm 20Q^. 
Although the fre q uency of PISNe is debated (e.g. , 
Scannapieco et a l.l 120 031: iVenkatesan &: TruranI 120031 : 
Norman et al.l 12004: Ri cotti fc Ostrikeii I2004D . we have 
found that the overall dynamics of the SN remnant 
and the distribution of metals are largely indepen- 
dent of the progenitor and are governed mainly by 
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Fig. 14. — Hydrogen number density averaged along the line of sight in a slice of IO//1 kpc (comoving) around the x-y plane, overlaid with 
the distribution of all metal particles (bright orange), after 1, 10, 50, and 200 Myr. The extent of each particle is set by the SPH smoothing 
length, while the apparent mixing within neighboring minihalos is a projection effect, (a): After 1 Myr, the shock has left the host halo 
and enters the IGM, while the metals are still confined well within the virial radius, (b): The contours of the shock have become visible as 
the SN remnant plows into the IGM 10 Myr after the SN explosion. The interior bubble expands adiabatically into the cavities created by 
the shock and begins to lose its spherical symmetry, (c): After 50 Myr, the dense shell becomes distinguishable from the metal-enriched 
interior, which has become substantially asymmetric and expands into the voids around the y-z plane in the shape of an hourglass, (d): 
When the SN remnant finally stalls after 200 Myr, the dense shell and the metal-enriched, interior bubble have reached their maximum 
extent, with the former providing a natural confinement around the interior. 



the explosion energy. Our simulations would there- 
fore also approximately describe a hypernova explo- 
sion, where a rotating, massive star undergoes core 
collap se (e.g., lUmeda fc Nomotol 120021 : iTominaga et al.l 
[2001 . Such scenarios might better explain the pecu- 
liar yields found in extremely metal-poor stars in the 
Galactic halo (e.g., iChristlieb eFall l2002t iFrebel et all 
[2005h . assuming they formed out of gas enriched by 
a previous genera tion of stars (e.g.. lUmeda fc Nomotol 
2003, 2005; TumU nson et al.ll2004t i lwamoto e t al. 2005; 
KarlssonI I200Q : ITumlinsonI I2006D . However, in light of 



the tentative iden tification of SN 2006gy as a PISN 
(jSmith et al.l 12007b . we find that the plausibility of the 
PISN scenario is enhanced and that it provides a viable 
possibility for the ultimate fate of a very massive, metal- 
free star. 

When and where do the first Pop II stars form? Based 
on the progression of structure formation, we expect that 
the first low-mass stars will form at the centers of the first 
galaxies, where primordial material streams in along fila- 
ments and mixes with the metal-enriched gas of the host 
halo. On the other hand, metals expelled into the voids 
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arc not so readily available, as the associated gas resides 
at too low densities and high temperatures. Because of 
this diversity, secondary star formation in the virializa- 
tion of the first galaxies will likely be a highly compli- 
cated and multi-faceted process. To understand the rel- 
evant mechanisms in detail, one must perform numerical 
simulations that include metal cooling as well as an ef- 
ficient model for their mixing. Furthermore, one must 
take into account the effects of multiple star-forming re- 
gions, including the expansion of additional H II regions 
and SN remnants. In future work, we hope to shed some 
light on these issues by performing detailed numerical 
simulations incorporating all the necessary physics and 



evolving them to sufficiently late times. Such investiga- 
tions are key in arriving at detailed predictions for the 
properties of the first galaxies to be observed with JWST. 
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